Non-Markovian Dynamics and Entanglement of Two-level Atoms in a Common Field 



C. H. Fleming* and N. I. Cummingsf 
Joint Quantum Institute, University of Maryland, College Park, Maryland 20742-4111, USA 

Charis Anastopoulos^ 
Department of Physics, University of Patras, 26500 Patras, Greece 

B. L. Hu § 

Joint Quantum Institute and Maryland Center for Fundamental Physics, 
University of Maryland, College Park, Maryland 20742-4111, USA 
(Dated: Jan 11, 2011) 

We derive the stochastic equations and consider the non-Markovian dynamics of a system of 
multiple two-level atoms in a common quantum field. We make only the dipole approximation for 
the atoms and assume weak atom-field interactions. From these assumptions we use a combination 
of non-secular open- and closed-system perturbation theory, and we abstain from any additional 
approximation schemes. These more accurate solutions are necessary to explore several regimes: 
in particular, near-resonance dynamics and low-temperature behavior. In detuned atomic systems, 
small variations in the system energy levels engender timescales which, in general, cannot be safely 
ignored, as would be the case in the rotating- wave approximation (RWA). More problematic are the 
second-order solutions, which, as has been recently pointed out [ ], cannot be accurately calculated 
using any second-order perturbative master equation, whether RWA, Born-Markov, Redfield, etc.. 
This latter problem, which applies to all perturbative open-system master equations, has a profound 
effect upon calculation of entanglement at low temperatures. We find that even at zero temperature 
all initial states will undergo finite-time disentanglement (sometimes termed "sudden death"), in 
contrast to previous work. We also use our solution, without invoking RWA, to characterize the 
necessary conditions for Dickie subradiance at finite temperature. We find that the subradiant states 
fall into two categories at finite temperature: one that is temperature independent and one that 
acquires temperature dependence. With the RWA there is no temperature dependence in any case. 



I. INTRODUCTION 

A. Motivation 

Atomic systems constitute an important setting for the 
investigation of quantum decoherence and entanglement 
phenomena essential for quantum information processing 
considerations [2—5]. The physical principles underlying 
these systems are quite well understood, and they can 
be controlled and measured with great precision. One 
aspect of quantum entanglement dynamics that has re- 
ceived significant attention is the phenomenon of entan- 
glement sudden death, or finite-time disentanglement, 
while energy and local coherences only decay away expo- 
nentially in time [6-8] . A common setting for theoretical 
discussion of this phenomenon is atomic systems inter- 
acting with the electromagnetic field [ , 9-13], serving 
as an environment in the quantum open system (QOS) 
perspective. 

When considering atomic dynamics with an open- 
system approach, where there is a continuum of field 
modes that are treated as a reservoir, the Born-Markov 
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approximation (BMA) [14-17] is commonly invoked. The 
BMA is sufficient to calculate the second-order (in the in- 
teraction Hamiltonian) dynamics of quantum open sys- 
tems, even when far outside of the Markovian regime 
[17], such as low-temperature and 1// noise. However, 
it has been recently discovered that the second-order dy- 
namics (more specifically, the master equation) of quan- 
tum open systems is not sufficient to fully generate the 
second-order solutions, including detailed state informa- 
tion such as entanglement [1]. The errors involved are 
particularly exacerbated at low temperature where the 
noise correlations of the environment exhibit long-ranged 
correlations, whereas there are no such problems in the 
Markovian regime wherein the BMA is exact. Therefore, 
when calculating entanglement dynamics with a second- 
order master equation, careful attention must be payed 
to whether a positive value of the entanglement measure 
exceeds the inherent second-order errors. 

Much of the theoretical work on atom-field systems 
is also derived using the rotating-wave approximation 
(RWA) [14, 15, 18, 19], which eases calculation by plac- 
ing the master equation into a Lindblad form [20-31] 
and therefore ensuring completely-positive evolution. In 
addition to weak coupling the (post-trace) RWA also re- 
quires that there be no nearly-resonant energy levels 1 . 



Near-resonant terms can be preserved in implementing the RWA, 
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As most previous work on multipartite open systems has 
been performed at resonance, it is worth asking how 
strongly such results are dependent upon this assump- 
tion. 

Until now, a systematic treatment of multiple two-level 
atoms in a common quantum field has yet to be carried 
out in a manner which can fully and accurately predict 
entanglement evolution and address critical issues such as 
sudden death of entanglement. There are two important 
reasons for this. (1) Calculations have usually involved 
perturbative master equations, either explicitly or by in- 
voking the BMA or RWA which both implicitly assume 
a perturbative coupling to the environment 2 . However, 
second-order master equations fundamentally cannot rule 
out second-order errors in their entanglement monotones. 
This is due to the little-known fact that second-order 
(non-Markovian) master equations are not generally ca- 
pable of providing full second-order solutions except at 
early times [1]. Such inescapable errors can potentially 
lead to entirely different and erroneous qualitative fea- 
tures of entanglement dynamics and sudden death. (2) 
Use of the RWA also does not allow for consideration 
of near resonance (as additional near-stationary terms 
are needed in the Dirac picture) The existence of a sub- 
radiant dark state generically requires the resonance con- 
dition, but determining how critical this is requires some 
analysis of the near-resonance regime. The aim of the 
paper is to lead towards a complete and systematic treat- 
ment of multiple two-level atoms in a common quantum 
field. To this end, the issues (1) and (2) have to be ad- 
dressed and this is what we do in this paper. 



B. New Methods and Results 

In this work we employ the mundane second-order 
master equation, which we describe in Sec. Ill, but only 
to explore the open-system dynamics and zeroth-order 
state information (e.g. dark and bright state forma- 
tion). At second order the perturbative master equation 
is consistent with the well-known Born-Markov and Red- 
field master equations, however it does not require any 
Markovian approximation to be derived [17, 33-35]. We 
generally forsake the RWA so that we can explore near- 
resonance, however this is otherwise irrelevant in deter- 
mining this timescale information[19]. The new method 
we apply in utilizing the well-known second-order master 
equation are the canonical perturbative solutions detailed 
succinctly in [1] and more thoroughly in [17], which we 
also briefly outline in this paper. 



but this will then lead to a master equation not of Lindblad-form 

as in [32]. 

2 We would note that a non-perturbative (or higher-order) treat- 
ment of the model with an RWA system-environment interaction 
is still qualitatively interesting. 



To calculate the late-time entanglement of the open- 
system, as is vitally important to consider questions per- 
taining to death of entanglement, we require second-order 
state information which no second-order master equa- 
tion can provide [.!]. In principle one could determine 
this information from the fourth-order master equation, 
which naturally excludes the BMA and RWA, however we 
choose to more simply calculate the second-order asymp- 
totic state directly from (canonical) closed system + en- 
vironment perturbation theory in Sec. V. Note, there- 
fore, that this second-order state information is necessar- 
ily generated by non-Markovian dynamics, even though 
we avoid such an open-system calculation. We are able to 
calculate the second-order asymptotic state of the open 
system at zero temperature and at high temperatures. 
This is a new method that we have developed specifi- 
cally to overcome shortcomings of the perturbative mas- 
ter equation discovered in [1], and more attention will be 
given to it in future letters. 

In this way we are able to show that the two atoms in a 
single field are not asymptotically entangled, even when 
near resonance and very close together — which is the 
criterion for a dark state. This asymptotic behavior turns 
out to be rather opposite to that of two oscillators in a 
field, which can be asymptotically entangled [36]. In fact, 
we find that even at zero temperature the entanglement 
of any pair of atoms will always undergo sudden death, 
regardless of the initial state. In addition to this we also 
find that the prediction of revival of entanglement in [11] 
can be unreliable in some cases (due to similarly-sized 
errors arising from the use of the second-order master 
equation), unless the decay rate is sufficiently small. 

Furthermore we explore the existence of sub and super- 
radiant (dark and bright) Dickie states at finite temper- 
ature for many atoms and with detuning for two atoms 
in Sec. IV B. To achieve dark and bright states one re- 
quires proximity better than the resonant wavelength, as 
is known, and tuning better than the ordinary dissipa- 
tion rate, as we find. The effect of finite temperature 
is relatively more interesting, especially when one does 
not employ the RWA-interaction approximation in the 
Hamiltonian. There appear to be a subset of dark states, 
which we designate improper dark states, which become 
temperature dependent. 

We now briefly discuss the Markov property and its 
implications upon open-system dynamics. In Sec. II 
we present the microscopic model which generates our 
stochastic equations of motion. In Sec. Ill we briefly 
describe the perturbative second-order master equation. 
Then in Sec. IV we discuss the corresponding master- 
equation solutions, the resulting dynamics, and the ac- 
curacy of these solutions. Finally in Sec. V, we discuss 
the alternative calculation of the second-order asymp- 
totic state, as well as the implications upon entanglement 
dynamics. 
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C. Non-Markovian Dynamics and 
Completely-Positive Evolution 

To clarify our language, we use the terminology Marko- 
vian and non-Markovian in the classical sense. The en- 
vironment and resulting dynamics of an open system 
are said to be non-Markovian if the underlying stochas- 
tic processes of the environment exhibit non-singular 
time correlations or non- vanishing correlation timescales. 
In this sense, Markovian dynamics correspond to white 
noise whereas non-Markovian dynamics correspond to 
colored noise. For a relativistic field, the noise is col- 
ored by both low-temperature (quantum) fluctuations 
and the (relativistic) dispersion relation (see Sec. II B). 
This is the more physical notion of the Markov property, 
as given non-Markovian dynamics in this sense (1) the 
quantum regression theorem does not apply [37], (2) one 
cannot add Hamiltonian terms to the master equation 
post derivation (e.g. see [38]: Sec. 7.1), (3) one cannot 
add dissipative terms to the master equation post deriva- 
tion (e.g. see [32]). 

The classical usage of the terminology Markovian and 
non-Markovian is slightly different when applied to the 
master equation itself (and not its resulting dynamics). 
A master equation is said to have a non-Markovian repre- 
sentation when placed into integro-differential form, e.g. 



p(t) = f drK{t-T)p{j), 
Jo 



(1) 



whereas the master equation is said to have a Markovian 
representation when placed into a time-local form, e.g. 



p(t) = C(t) p(t) 



(2) 



Obviously this distinction is superficial here as one can 
trivially relate these two master equations 



c(t) = g(t)g(t)-\ 



Q{s) 



s — K.(s) 



(3) 
(4) 



in terms of the Laplace transform, defined /(s) = 
^dte~ st f(t). While Markovian dyn amies imply a 
Markovian representation, the inverse is not always true, 
e.g. see the Hu-Paz-Zhang master equation [38, 39]. 
Time-dependence in C(t) is also not specific to non- 
Markovian dynamics, as one can couple a time-dependent 
Hamiltonian system to a Markovian environment. 

These distinctions are further complicated for Lindblad 
master equations obtained after the post-trace RWA [19]. 
Given their mathematical structure, such master equa- 
tions can be unraveled with Markovian processes, even 
if the microscopic environment from which they were de- 
rived was not Markovian (or even close to being Marko- 
vian). This Markovian theory is essentially an effective 
theory, which is sufficient for modeling certain aspects of 
the perturbative dynamics. 

One might question how positivity can be maintained 
with a master equation which is not of Lindblad form. 



In general, for non-Markovian dynamics, master equa- 
tions need not be of Linblad form to ensure completely- 
positive evolution. Indeed, there are exact master equa- 
tions known which do not have Lindblad form and yet 
have a time-local or Markovian representation, e.g. see 
the Hu-Paz-Zhang master equation [38, 39]. The positiv- 
ity, or lack thereof, generated by a perturbative master 
equation not of Lindblad form can be understood as fol- 
lows. If we can produce solutions p(t) which are good to 
some order 0(g 2n ) in the system-environment coupling, 
parameterized by <?, and the higher-order errors are non- 
secular in time, then positivity is only ensured to that 
order, or more specifically 



i^p{t)i) > + O(g 



2n+2\ 



(5) 



for all Hilbert-space vectors ip. Manipulating the mas- 
ter equation into a Lindblad form will produce solutions 
which satisfy Eq. (5) exactly, however they are not nec- 
essarily more accurate in any other sense. The pseudo- 
Lindblad master equation will preserve the trace and Her- 
miticity of the density matrix exactly, and the master 
equation will belong to the appropriate class of algebraic 
generators (being the difference of two Lindblad genera- 
tors) . 



II. THE MICROSCOPIC MODEL 

A. Interacting Hamiltonians 

We wish to investigate the properties of multiple atoms 
interacting with a common electromagnetic field in free 
space, which serves as the environment in the open quan- 
tum system description. We will consider the atoms to be 
held stationary at fixed positions r„, and only consider 
the dynamics of the electrons in orbit. As we consider 
neutral atoms far apart, we will neglect the electrostatic 
interactions between the atoms. We will use the two- 
level approximation to describe the atoms, so that they 
are an array of, otherwise non-interacting, qubits. Only 
the lowest-level excitations of the Dirac field are consid- 
ered (see for instance, [ ■ ], App. A), so that we have 

H sys+ fi id = H sys + Hi n t + Hfioid , (6) 

Hsys = ^ °+„ CT - n j (7) 



H; 



(8) 



where er denote the Pauli matrices and d corresponds to 
the lowest-energy levels of the electron's dipole matrix 
[18, 40, 41]. The field Hamiltonian and vector potential 
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are given by 



H fi eid = /// d 3 k oj k ak 6i a ktCh , 



A(x) 



1 



(2tt) 3 / 2 



d 3 fc^A k , e)c (x). 



(9) 
(10) 



{■ 



+ ik-x 



ak, £ 



(11) 

with k — ujk/c and where Ek denote the polarization vec- 
tors perpendicular to k [42, 43]. We have assumed the 
atoms to be relatively stationary and that the atomic 
transition in each atom will produce linearly polarized 
photons (i.e., both ground and excited state are eigen- 
states of some component of angular momentum with 
the same eigenvalue). 

Even though we assume the atoms to be ideal dipoles, 
we do neglect the 1/r 3 electrostatic interactions between 
the electronic dipole moments, and treat r as a classi- 
cal parameter rather than a quantum degree of freedom. 
Thus we necessarily assume a sufficient degree of separa- 
tion between all atoms. It will be seen that magnetostatic 
interactions between the atomic dipole moments will be 
introduced via backreaction from the field, but these will 
be terms will be 1 jr and thus can play a role even in the 
regime where the electrostatic interactions are negligible. 



B. Environment Correlations 

Any perturbative open-system analysis of this prob- 
lem, whether the master equation we employ here or 
Langevin equation, will entail calculation of the corre- 
lation function for the field we couple to 



a nm (t,r) = (d£ A(r„;f) A(r m ; r) T d m ). 



field 



(12) 



where A(r; t) denotes the time-dependent vector poten- 
tial A(r) in the interaction or Dirac picture. All per- 
turbative master-equation coefficients in Sec. Ill will be 
generated by these free-field correlations. In addition to 
being Hermitian and positive definite, for our considera- 
tion of stationary atoms, the correlation function is also 
stationary: a nm (t,T) = a nm (t — r). For simplicity we 
will further assume that all system dipole moments are 
of the same magnitude in this work, though such differ- 
ences can be easily accounted for. 

The fluctuation-dissipation relation (FDR) [44] allows 
us to express the environmental correlations in terms of 
the damping kernel as 



a(oj) = 7(w) 



sinh(^) 
27(w)w[n(M,T)+0(-w)] , 



here in the Fourier domain, defined 



dti 



fit) 



(13) 
(14) 

(15) 



and where n(oj, T) is the thermal average photon number 
in a mode of frequency oj. 

The damping kernel is formed by commutators of the 
vector potential and therefore in these models it is inde- 
pendent of the state of the environment (temperature in- 
dependent) and also the same whether in the classical or 
quantum regime. Outside of the high-temperature (semi- 
classical) regime where a(u>) = 2Tj(u), the damping 
kernel is the simpler object to inspect, and in the free- 
particle Hamiltonian corresponding to (6) it directly pro- 
duces the Abraham-Lorentz force [45]. At least petur- 
batively, the damping kernel is only responsible for irre- 
versible dissipation of energy, whereas the full correlation 
also contains the influence of noise [44] . 



1. Scalar-Field Correlations 

First let us consider the simpler case of a scalar field, 
wherein one essentially neglects the polarization vectors 
in A and the direction of atomic dipole moments. In this 
case we calculate the associated damping kernel to be 



7nm(w) = 70 sinc(r„ m w) , 



(16) 



where r nm = r„ — r m is the difference vector for the po- 
sitions of the two atoms and 70 = TVmO^-O is the self- 
damping coefficient. Restoring factors of c, the sine func- 
tion appears to act as a relativistic regulator for oj <C c/r. 
Therefore, in the nonrelativistic regime one only has lo- 
cal damping for charges at precisely the same location. 
In fact, this must always be the case for linear coupling 
to a relativistic field with a position-like field operator 
such as A. The damping kernel is deterministic and the 
same whether in the classical or quantum regimes. Thus 
it is necessarily constrained to the light cone. 

Note for instance the temporal representation of the 
scalar-field damping kernel 



Jnm(t) = y6r„ m (i) 

c /.a _ 0(r-\t\) 



(17) 
(18) 



where 8 is the Heaviside step function. This kernel 
strictly adheres to the light cone. 



2. Electromagnetic- Field Correlations 

The electromagnetic damping kernel is slightly more 
complicated as it contains orientation dependence from 
the polarization vectors. We calculate this kernel to be 



7r. 



% {oj) = 70 d„ [FSi(r nm oj) + FS (r nm oj) r nm ij im } d m 



7o 



1 

6"7T£oC 3 



(19) 
(20) 
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FIG. 1: Comparison of sine (bold), FSi, and FS 
(dashed). Sine and FSi are extremely qualitatively 
similar, both being unity at zero whereas FSo vanishes 
at zero. 




FIG. 2: Comparison of the same functions in Fig. 1, but 
in the time domain: the rectilinear distribution (bold), 
FSi, and FSo (dashed). All of these functions are 
constrained to the light cone. 



in the Fourier domain. Instead of a sine functions, we 
have the related entire functions 



3 (z 2 - l)sin(z) + zcos(z) 
tbi(z) = +- -= 



FSq(z) 



3 (z 2 — 3) sin(z) + 3 z cos(z) 
2 ^ 



(21) 
(22) 



In Fig. 1 we compare these functions. In any case, the 
cross correlations are very nonlocal and are maximized 
when the atoms are close. Whereas when the atoms are 
very far apart, the cross correlations always vanish and 
thus all noise can be treated independently. One can 
see that the scalar-field correlations are very similar to 
that of the electromagnetic field when d„ || d m _L r nm . 
As we will wish to maximize cross correlations, we will 
primarily work with the scalar-field correlations, which 
are very similar to the problem of parallel dipoles in the 
electromagnetic field. 



III. SECOND-ORDER MASTER EQUATION 

A. Master equation and coefficients 

The second-order master equation for the reduced den- 
sity matrix of the dipoles can be expressed [J 7] 



p = £p, 
c = c + c 2 



(23) 
(24) 



in terms of the zeroth and second-order Liouville opera- 
tors 

£oP=[-*H,p], (25) 

£2 P = ^ l^™ ' P ("A- m nO CT Xm ) f - (An m O <T Xm ) p\ . 

nm 

(26) 

We present the stationary limit of this master equa- 
tion for the scalar- field environment, by way of (27), 
(29), (13), (16), (31) (finite temperature) and (33) 
(low-temperature expansion). Some aspects of the 
electromagnetic-field environment and full-time theory 
are also discussed herein. We choose to present more re- 
sults for the scalar-field environment merely as a means 
of presenting simpler expressions which retain the quali- 
tative features we are interested in. The electromagnetic 
theory can be resolved without a considerable amount of 
additional difficulty, but the resulting expressions will be 
several times more lengthy. 

The second-order operators in master equation (26) are 
most easily represented by the ladder operators as 

(27) 



1 r i 



(28) 



and the second-order coefficients being generated from 
the field correlations as 



A nm (u>) = ~a nm {u) - iV 



* Or, 



(29) 



here in the late-time limit (as compared to system and 
cutoff timescales), where V denotes the Cauchy princi- 
pal value. Higher-order master equation coefficients will 
entail convolutions over several copies of the field cor- 
relations combined with several products of the system 
coupling operator. 

The first portion of the second-order coefficient, or Her- 
mitian part (here real), is immediately given by Eq. (13). 
Whereas the second term, or anti-Hermitian part (here 
imaginary) , must be evaluated via the convolution 



Im[A lm M] = - 



1 

2^ 



deV 



1 



(e) , (30) 



G 



and together they form a causal response function. These 
are the coefficients which often require regularization and 
renormalization. For now let us simply evaluate the bare 
coefficients for non-vanishing r. For finite temperatures, 
the coefficients exactly evaluate to 



1m[Anm(u)] 



7o 1 



fin 



coth( — ) - I 
2TJ 



cos(r m „cj) 



in terms of the Lerch <&i function 



-A- A 



fc=l 



k + z 



(31) 



(32) 



This functional representation is exact, though best for 
positive temperature. Conversely, one also has the low- 
temperature expansion 



x r a ( si 7o sign(w) ^ c 



(33) 



7o 1 



sin(r rim w) ci(|r„ ro w|) - cos(r nm w) si(r nm w)] 



'nm ?f 

in terms of the summand 

Ei[(+fc/3 + ir„ m )|w|] Ei[(-fc/3 + ir„ m )|w|] - ztt 



,(+kfi+ir nm )\u\ 



o(-k0+ir nm )\u]\ 



and where the trigonometric integrals are defined 

,sm(z') 



si(z) 
ci(z) 

Ei(z) 



dz 

o 

dz 



z' 

,cos(z') 



dz'V 



(34) 

(35) 
(36) 

(37) 



however, for positive temperatures this expansion is not 
well behaved for small energy differences. For zero tem- 
perature, the exact relation (the second line in (33)) 
is well behaved and matches perfectly to the zero- 
temperature limit of Eq. (3f ). 



B. Unitary Dynamics and Renormalization 

First we would like to focus upon the unitary gener- 
ators that have been induced by the environment. If 
we express the master equation in the standard pseudo- 
Lindblad form 

£ 2 {p} = hU,p]+D{p}, (38) 
then the unitary generator U is given by 

U = \ a *n (Aim<> <T X ) - (Am* Vx„y CTxA ■ 

2,1 * — ' 

nm 

(39) 



One must be careful not to identify all of U with 
environmentally-induced forces, such as dissipation and 
renormalization. The delineation between homogeneous 
backreaction and noise is not as clear in the master equa- 
tion as it is in the Langevin equation (the Heisenberg 
equations of motion for the system operators as driven by 
the field). The Langevin equation instructs us that only 
the self interactions generated by Im[^4„„(0)] = — 7„ n (0) 
should be renormalized [17, 45]. 

From the master-equation perspective this prescription 
can also be motivated, though it is more difficult to im- 
mediately see. The self terms U„„ happen to be the most 
divergent. These terms resolve to be 

U n „ = lm[An n (+Q, n )] <J- n er +n + ImL4„ n (-fi n )] o+„ °"-„ 

(40) 

The most divergent part of this expression is not present 
in the energy-level separation between the ground and 
excited states, but only in the average or absolute of the 
energy levels. Therefore if we choose the renormalization 
suggested by the Langevin equation 

U rcn = Vxn Im[Am(0)] cr Xn , (41) 

n 

then the divergent contribution from 



Im[Am(0)] = - 



7o 

2 T r ,r 



(42) 



will indeed be removed. 

It is also physically instructive to inspect the corre- 
sponding cross terms from Im[.A„ m (0)], as the Langevin 
equation reveals these terms to be non-dissipative back- 
reaction mediated by the field (with Im[^4„„(0)] being 
their corresponding divergent self interaction). For the 
electromagnetic field these terms can be assigned the in- 
teraction potential 



, r M0 V"^ ,Tl + r «»Ti ] 

V 2LA = 2^ d « r 



(43) 



which is qualitatively approximated by the scalar field 
interaction for d„ || d m _L r nm . This result for the 
two-level atoms should be compared to the well-known 
magnetostatic potential of the (classical) Darwin Hamil- 
tonian for point particles 



V Da 



Mo 
8tt 



E 



Pn 



1 



(44) 



and so these two-level atoms experience magnetostatic 
attraction and repulsion depending upon their relative 
momentum states, just as moving point charges do. 



1. Regularization 

Although we have renormalized away the 1 / r nn diver- 
gences present in Im[^4„„(0)], there remain logarithmic 
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r nn divergences in ImL4„„(f2„)] which cannot be renor- 
malized, as they extend into the dissipative portion of 
the master equation. These divergences must be made 
finite with regularization. 

Note that sinc(w/A) is a high frequency regulator: 
sinc(z) : [0, oo) — > [1,0) sufficiently fast for all of our 
integrals to converge. Therefore we don't need to con- 
sider any additional regularization in our damping kernel 
if we do not evaluate sinc(rw) for vanishing r. Instead 
of allowing r to vanish for self-correlations, we impose a 
high frequency cutoff r = A -1 , perhaps motivated by 
the non-vanishing physical size of the dipole. The more 
common alternative is to introduce cutoff regularization 
directly into the field coupling in Eq. (10), often by 
treating the coupling strength as a form factor with some 
gradual k-dependence. Different choices of cutoff regu- 
lators will yield the same results to highest order in A, 
and the theory should be somewhat insensitive to these 
details in the end. 



C. Full-time non-Markovian dynamics 

For the full-time evolution of initially uncorrelated 
states, one must apply the full-time coefficients (again 
we work from the notation of Ref. [17]) 

Anm(w,t)= [ dr e~ lulT a nm {f) , (45) 
Jo 

which exhibits causal behavior in as much as the field 
correlations are causal. In the semiclassical regime the 
field correlations are perfectly causal as a. = 2T~f, and 
the damping kernel 7 is classical and perfectly causal 
(see Fig. 2). However, in the quantum regime low- 
temperature fluctuations are allowed to leak outside of 
the light cone. 

If we start with factorized initial conditions for our 
atoms and field, then classically they would have no 
knowledge of each other's existence until the mediating 
photons could reach them at t = r/c. Thus the dynam- 
ics would be perfectly factorized until this time (a highly 
non-Markovian behavior). In the quantum regime, the 
story is slightly different as unavoidable quantum fluctu- 
ations do leak from one atom to the other, in a manner 
which only roughly adheres to the light cone. We plot the 
resulting dynamics of the coefficients in Fig. 3, where one 
can see the master-equation coefficients jolt at precisely 
t = r/c. Prior to the jolt, the master equation coefficients 
(for finite r) are roughly zero and the dynamics factor- 
ized; whereas after the jolt, the coefficients are roughly 
their asymptotic value and dynamics interrelated. 

Shortly after t — r/c, one has what is essentially a joint 
state of the system which is properly correlated with the 
field. Properly-correlated initial states can be more accu- 
rately modeled by a more laboratory-appropriate prepa- 
ration scheme wherein a desired state of the open sys- 
tem is properly prepared from a global equilibrium state, 



yet without significantly disturbing the vital system- 
environment equilibrium correlations [ ]. However in 
the weak-coupling regime it is sufficient to consider the 
asymptotic dynamics. Therefore from hereon we will con- 
sider jC(oo) as approximately corresponding to the open- 
system dynamics of properly-correlated initial states, and 
not the highly non-Markovian dynamics of factorized ini- 
tial states. 



IV. SECOND-ORDER DYNAMICS 
A. Canonical perturbation theory 

The open-system dynamics of properly-correlated ini- 
tial states are approximately described by the time- 
independent Liouvillian £(00), which we will now write 
simply as C The time evolution is then approximately 
e tc , and it can be computed (analogously to the time- 
independent Schrodinger equation) simply from the so- 
lutions of the eigen- value problem 

£{o} = /o, (46) 

where / is an eigen-frequency and o a right eigen- 
operator (super- vector) 3 . In principle this can be per- 
formed numerically with the super-matrix operators, but 
to avoid more costly numerics we have resorted to a care- 
ful application of canonical perturbation theory, as can 
found in Ref. [17]. Because the master equation itself is 
perturbative, there is no loss in accuracy by finding the 
solutions perturbatively. 

It is of paramount importance to note that second- 
order master equations can only determine all eigen- 
values / to second order [1]. This includes the coherent 
oscillation, emission, and decoherence rates. However, 
entanglement and other state information is obtained 
from the eigen-operators o, which second-order master 
equations can generally only determine to zeroth order 
[1]. Therefore in this section we only address dynamical 
questions answered by /, and in Sec. V we will calculate 
and inspect p{oo) to second order, by using an alternative 
mathematical formalism which has similar ingredients. 

At zeroth-order our eigen-value problem corresponds 
to the energy-level differences and outer-products of en- 
ergy states 

£o{|wi)(w i |} = -iLOij \oJi)(oj 3 \ , (47) 

where u>ij = — Lij . The environment induces frequency 
shifts (including decay) and basis corrections such that 
the eigen-operators are no longer dyadic in any basis of 
Hilbert-space vectors. Some degree of degeneracy is also 
inescapable as to a = uijj = 0. 



3 These eigen-operators are often referred to as the damping basis 
of the master equation [4 7]. 
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FIG. 3: Re[.A„ m (— f2 m ; t)] (left) and Im[./l nm (— f2 m ; i)] (right) for a zero temperature reservoir at r nm — 10/J7 m . The 

bold line denotes the asymptotic coefficients. 



As our system coupling is non-stationary, with no addi- 
tional degeneracies the cross-coupling will have no effect 
upon the second-order frequencies of the perturbed off- 
diagonal operators, and the /y corresponding to |u>j)(wj| 
for i 7^ j are given by 



fij = -lUli 



£ 2 {|w i )(w i |} \tjj 



(48) 



which reference no cross-correlations. Second-order cor- 
rections to the eigen-operators o (and thus states) can 
then be found by perturbative consistency with the mas- 
ter equation. Dynamics of the diagonal operators and 
any other degenerate (and near-degenerate) subspaces 
must be treated much more carefully with degenerate 
perturbation theory. For the energy states, their second- 
order dynamics are encapsulated by a Pauli master equa- 
tion. This gives rise to their second-order relaxation rates 
and zeroth-order eigen-operators. Due to inherent degen- 
eracy, uJa = ojjj = and any resonant frequencies, their 
second-order operator perturbations require the fourth- 
order Pauli master equation [1, 17]. To summarize, in 
general the matrix elements of the solution p(t) expressed 
in the (free) energy basis will be accurate to 0(7) off 
the diagonal but only to 0(7°) on the diagonal (though 
timescales are known to 0(7)). This inaccuracy in the 
diagonals is an inherent limitation of any perturbative 
master equation, including those derived under the RWA 
or the BMA [ ]. With the RWA, however, all matrix 
elements are only good to O(-f ) 4 . 

In Figs. 4-5 we plot all relaxation rates associated with 
the two-atom system as a function of proximity, where 
7 is specifically the decoherence rate of a single isolated 
atom. These are quantities that can be calculated from 
the RWA-Lindblad equation, and our results are consis- 
tent with those reported in [11]. For large separation the 



4 When looking only at observables time-averaged over many sys- 
tem periods 2iz/Q some of these additional discrepancies gener- 
ated by the RWA can be greatly reduced. 
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|0.0> 



FIG. 4: Decay rates of the (zeroth-order) stationary 
operators for two resonant atoms in a zero-temperature 
environment at varying separation distance. The legend 
indicates the pure states they approximately correspond 
to in the order they occur at the vertical axis. 




FIG. 5: Decoherence rates of the (zeroth-order) 
non-stationary operators for two resonant atoms in a 
zero-temperature environment at varying separation 
distance. The legend indicates the matrix elements they 
correspond to in the order they occur at the vertical 
axis. 



decay rates for |%) = (|0,1) ± 1 1,0) ) js/2 are 1 + 1 times 
7 (which would be A7 for N atoms), as the noise pro- 
cesses are independent and the decay rates arc additive. 
Whereas at proximity they become and A 2 times 7 
for |\&_) and \*f?+) respectively, as the noise processes are 
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FIG. 6: Decoherence rates of the (zeroth-order) 
non-stationary operators for two atoms in a 
zero-temperature environment at varying detuning and 
vanishing separation, 7*12 <C f2i,f22, with 7 = (fi)/100. 

The legend indicates the matrix elements they 
approximately correspond to for small detunings (in the 
order they occur at the vertical axis) as to compare 
with Fig. 5. 



maximally correlated and display destructive and con- 
structive interference. 

In Fig. 6 we plot all non-stationary decoherence rates 
associated with the two-atom system as a function of de- 
tuning. These are new results that cannot be determined 
from an RWA-Lindblad equation, as for slight detuning 
there are near-stationary contributions to the dynamics 
(in the interaction picture) which cannot be discarded. 
To achieve a dark state, the tuning of the two atoms 
must be much better than the dissipation, SQ -C 7, which 
counter-intuitively implies that weak-dissipation is not 
always desirable to preserve coherence. However, this 
condition makes more sense if thought of in another way: 
The dark state arises from the destructive interference of 
the emission from the two atoms. If the emission from 
each atom is characterized by center frequency O n and an 
emission line width 7, then the condition 50 T simply 
specifies that the emission lines of the atoms must over- 
lap enough that their emissions are not distinguishable 
from one another. This allows the required destructive 
interference. 



1. The Atomic Seesaw 

One behavior which is qualitatively different from the 
closed-system evolution is the damped oscillations be- 
tween the singly-excited states. More specifically for any 
initial state of the form 



Po 



10,1} 
|1,0> 



a + S 
-ib 



+ib 
a — S 



(OTI 
(1,01 



(49) 



with all positive coefficients, then in addition to the Bell- 
state decay one will also have damped oscillations of the 



pcoB(At) + 6sin(/it)] e-^ 1 (|0,1)(0,1| 
-% [6cos(/ii) - *sin(/it)] e-^ 1 (|0,1)(1,0| 



|i,o)(i,o|) 
|i,o)(o,i|) 

(50) 



which can oscillate from one excited state to the other 
excited state with the frequency 



fi 



cos(fir) 
2 70 — 



(51) 



for all temperatures. As this oscillation rate becomes ar- 
bitrarily large with proximity, it should be physically ob- 
servable in some systems. Though of course, this model 
does break down for sufficiently small separations. 



B. The dark states 

First we will review the normal conditions for which 
dark and bright (sub and super radiant) states emerge 
for two atoms in this model. Then in Sec. IV B 1 we will 
discuss how at finite temperature, and without a RWA- 
interaction Hamiltonian, two classes of (Dickie) dark 
states emerge for TV atoms, which we designate proper 
and improper. 

All stationary (and thus decoherence-free) states p D of 
the open-system must satisfy the relation 



£p D = o, 



(52) 



and are thus right eigen-supervectors of the Liouvillian 
with eigenvalue 0. As the Liouvillian is not Hermitian, 
there is no trivial correspondence between the left and 
right eigen-supervectors. The super-adjoint of the mas- 
ter equation [16] time-evolves system observables and for 
closed systems can be contrasted 



C Q p = -i[H, p] , 
4s = +i[H,S]. 
if 



(53) 
(54) 



The left eigen-supervector S D corresponding to p D must 
therefore satisfy 



£ f S D = . 



(55) 



So for every stationary or decoherence-free state p D there 
is a symmetry operator Sd whose expectation value is a 
constant of the motion. The thermal state or reduced 
thermal state is such a state. In the limit of vanishing 
coupling strength, this state is the familiar Boltzmann 
thermal state. One can check that the symmetry op- 
erator in this case is proportional to the identity and 
corresponds to Tr[p] being a constant of the motion. 

For two resonant dipoles, with fl n = f2, there is an- 
other stationary state in the limit of vanishing separa- 
tion ri2 = r. Because of degeneracy, any superposition 
of states 



I*) =ai|l,0)+a 2 |0,l) 



(56) 
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is also an energy state and therefore annihilated by both 
Co and C Q . Further note that for vanishing separation, 
the field processes associated with A(r, t) become exactly 
correlated and identical. Their contributions to the in- 
teraction Hamiltonian can then be collected into 

H :i + H l2 = (<r Xl +rr X2 )d- A(r) = S x d- A(r) . (57) 

Next we note the equality 

E X |1,0) =E B |0,1) , (58) 

so that for the Bell states 

|%) = -L{|1,Q)±|0,1>}, (59) 

the noise adds destructively for |\IL) and constructively 
for 1 4'+). Therefore |^_) is a decoherence-free state (dark 
state) of the open system for vanishing separation and 
at resonance, regardless of coupling strength or temper- 
ature. And whereas |^_) appears dark (sub-radiant), 
appears bright (super-radiant). [Note that for anti- 
parallel dipoles, these roles will be reversed due to the 
anti-correlated noise.] 

In this particular case the left and right eigen- 
supervectors are equivalent, and so it is the dark-state 
component ( 1 I r _|p|4 r _) which is a constant of the mo- 
tion. However, unlike the thermal state, if the separa- 
tion is no longer vanishing then this is not some pertur- 
bative limit of a stationary state but of a very long-lived 
state. The final constant of motion, which we have val- 
idated by analyzing the eigen-system of C, corresponds 
to the coherence between the ground state and the dark 
state or (0,0 1 p |\&_). Using these constants of motion, for 
two very close dipoles in a zero-temperature environment 
with initial state p 0l the system will relax into the state 

Pl = (l-6)|0,OXO,0|+6|*_X*-| (60) 
+ c|0,0X*-|+c*|*_X0,0| , 

b= <tf-|po|*-> , (61) 

<0,0|p |*_) , (62) 

to zeroth order in the system-environment coupling, 
whereupon the system has bipartite entanglement b. 

While our (regularized) model is well behaved in the 
mathematical limit r — > 0, it is important to remember 
that physically the model is no longer valid for sufficiently 
small r. At small enough r other terms would come into 
play, including electrostatic interaction, and eventually 
the atoms would cease to even be distinct. We are assum- 
ing that this scale is much smaller than all other scales in 
our model (except perhaps the cutoff). This means that 
we can sensibly consider cases where r is small compared 
to the other parameters, but r cannot vanish completely. 

Since the coefficients of our master equation are contin- 
uous in r, it is useful to consider r = to understand the 
limiting behavior as r becomes small. The existence of 
the dark state we've discussed at r = means that this 
state will be almost completely dark when r is small; 



thus, any initial state p will first relax approximately 
into the state given in Eq. (60) within the ordinary relax- 
ation timescale 7, and then on a much longer relaxation 
timescale r, where roughly 1/r « 7(f2r) 2 for small r, the 
system will fully thermalize. However, this expression 
for the dark state is only to zeroth-order in the system- 
environment coupling. In order to understand the sub- 
sequent final state of decay one needs the second-order 
asymptotics that we discuss in Sec. V. 

Finally we would note that this "dark state" is a very 
general feature of resonant multipartite systems with 
similar linear couplings to a shared environment. One 
can rather easily work out that for a pair of resonant 
linear oscillators with these same noise correlations the 
sum mode is thermalized, and the difference mode is de- 
coherence free for vanishing separation. The separation 
dependence of the entanglement dynamics of two reso- 
nant oscillators was considered in Ref. [36] , while that of 
(effectively) two very close oscillators was considered in 
Ref. [48, 49]. 



1. Proper and improper dark states for N -atoms 

Here we specifically consider the finite-temperature 
theory of Dickie subradiance wherein we neglect all 
forces between the atoms, including the 1/r 3 electrostatic 
(dipole-dipole) and 1/r magnetostatic (orbital-orbital) 
interactions. It should be noted that for N > 2 the inclu- 
sion of atomic forces does make for significant differences 
in the theory which we do not consider here [50, 51]. 

The sub-radiant dark state achieves destructive inter- 
ference in the environmental noise (and thus littlc-to-no 
emission) while the bright state achieves constructive in- 
terference in the noise (and thus near-maximal emission). 
For the super-radiant bright state one essentially couples 
the system to N copies of the same field process associ- 
ated with A(r, t) and therefore the super-radiant emis- 
sion rate can be proportional to TV 2 . An N 2 dependence 
does appear the case as we demonstrate in Fig. 7. The 
emission rate is (perturbatively) determined by the noise 
correlation (the square of the noise process) . Both results 
differ having from N independent noise processes where 
one can simply add the N independent noise correlations 
which results in an emission rate at most proportional to 
N. Up to this point, the physics in this subsection is all 
well known. 

Following the previous approach in this section, we de- 
fine a proper dark state as an atomic state annihilated by 
Cq and Hi regardless of the state of the environment. Let 
us consider an assembly of N resonant dipoles at close 
proximity. We first note that the superposition 

I*) = ^2 a si,s 2 ,-,s N \si,s 2 ,--- ,s N ) , (63) 
£ " n =s 

of energy states with the same total excitement S is also 
an energy state and therefore annihilated by Cq. Defining 
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FIG. 7: Maximal (over states) second-order decay rates 
as a function of the number of atoms N, at zero 
temperature and in close proximity. The solid curve 
denotes the best quadratic fit and has a corresponding 
p- value of 2.4%, which is fairly significant in 
corroborating an A 2 dependence. 



the collective spin operator 



(64) 



such that the interaction Hamiltonian can be expressed 



H I = S x d-A(r): 



(65) 



a proper dark state must then satisfy Sa; I*-) = and 
will be decoherence free. For N — 2 this is the familiar 
Bell state that we've already labeled \^-)- 

In considering large iV the structure is essentially just 
what was studied by Dicke [52], so following that ap- 
proach we define collective y and z spin operators S y and 
S z as well as raising and lowering operators S + and 
analogously to Eq. (64), as well as S 2 = £ 2 + £ 2 + S 2 . 
And we can note that the free Hamiltonian for the atoms 
only differs from H z by a multiple of the identity, so 
all the eigenstates of that Hamiltonian are also eigen- 
states of S z . A basis for the Hilbert space of the sys- 
tem can be specified by the eigenstates of S 2 with eigen- 
values j(j + 1) and X z with eigenvalues m (though for 
N > 2 there will be degeneracy, so that additional quan- 
tum numbers are needed to identify a specific state). The 
dark state we seek must then satisfy S 2 I*-) = m \ ^ -) 
and H x \^f-) = 0. As the discussion in [52] implies, 
only states with j — and m — can satisfy these 
requirements simultaneously. Such states only occur 
when N is even, and that set of states has dimension 
Nl/ [{N/2 + l)!(iV/2)!]. These are also the dark states in 
the RWA, as they are in the null space of both S + and 
S_. For N = 4 these states take the form 

|*_) = ai(|0,0,l,l) + |1,1,0,0)) + 03(10,1,0,1) + |1,0,1,0» 
+ 03(10,1, 1,0) + |1,0,0,1)), (66) 

0=^>„, (67) 



where every pair in parenthesis is spin-flip symmetric. 
One can easily check that any such state is annihilated 
by £«. 

However, more generally we define an improper dark 
state as one only annihilated by C and not Hi (i.e., sta- 
tionary in the coarse-grained open-system dynamics but 
not in the full closed system dynamics), thus being de- 
pendent upon the state of the environment and even the 
coupling strength. In the simplest case we can consider 
the zero-temperature environment. For the second-order 
dynamics, upward transitions are automatically ruled out 
from the lack of thermal activation. Rather than investi- 
gating the master equation, for zero temperature we can 
then simply demand that the lowest-order decay transi- 
tions are vanishing, meaning that if \^^_) has total ex- 
citation 5", then (S'\ S a |**) = for all S' < S lesser 
and equally excited states. We can also state this in 
terms of the collective spin operators we have defined, 
by saying that we demand that | v I / -_) is an eigenstate 
of S z with eigenvalue m, and that all matrix elements 
onto states with lower m! values must vanish. Since 
Sjc = | (£+ + £_), we know that there will be non- 
vanishing matrix elements onto states with m! = m — 1 
unless 7ji — —j. So any state with m = —j is an im- 
proper dark state at zero temperature, and there are 
N\(2j + 1)/ [(A/2 + j + l)!(A/2 - j)!] such states [52]. 
Interestingly, for the RWA-interaction Hamiltonian such 
states (when combined with a vacuum field) are also sta- 
tionary states but of the closed-system dynamics. For 
N = 3 and at zero temperature, all such dark states can 
be expressed 

|¥_) = oi |1,0,0) + a 2 |0,1,0) + a 3 |0,0,1) , (68) 
= ^a n , (69) 

n 

for weak coupling to the field. Numerical investigation 
of the second-order master equation, and counting the 
number of null eigen-operators of C for different N and 
T, reveals that these dark states also exist for positive 
temperature, but they take on a different, temperature- 
dependent form. Though it is clear how the simple cri- 
teria leading up to resolution of the zero-temperature 
improper dark states fails to carry over into the finite- 
temperature regime, it is not clear how these finite- 
temperature improper dark states can be calculated with- 
out the aid of a finite-temperature master equation, nor 
what simpler properties they might exhibit. Moreover, 
they may only exist perturbatively in the second-order 
master equation: the fourth-order master equation might 
assign them fourth-order emission rates. 



V. SECOND-ORDER ASYMPTOTIC SOLUTION 

A. Canonical perturbation theory 

In this section we specifically address the asymptotic 
state p(oo) to second order. As we have previously men- 
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tioned, no second-order master equation can fully deter- 
mine all second-order state information, including mea- 
sures of entanglement. However, as we will point out, the 
second-order master equation is sufficient for determining 
the asymptotic state at high temperature. Moreover, at 
zero temperature one can apply canonical perturbation 
theory to the closed system + environment Schrddinger 
equation and bypass the master equation formalism en- 
tirely. 

To zeroth order in the system-environment interaction, 
the asymptotic steady state is Boltzmann, which can be 
expressed 



pt = n > 

Pr n = 



1 — tanh [ | a 



(70) 
(71) 



in terms of Pauli matrices. The asymptotic state of the 
second-order master equation is consistent with this re- 
sult and can additionally provide some of the second- 
order corrections 8p T via the constraint 



C {Sp T } + C 2 {p T } = . 



(72) 



These will specifically be the off-diagonal or non- 
stationary perturbations. In general, to find the second- 
order corrections to the diagonal elements of the den- 
sity matrix one needs to compute contributions from the 
fourth-order Liouvillian [1]. 

It has been shown [17, 53] that for non-vanishing in- 
teraction with the environment the off-diagonal elements 
of the asymptotic state match the reduced thermal state 



1 



Zc(j3) 



Tr, 



-0(H+H e +Hi) 



(73) 



where Zc{(3) is the partition function of the system 
and environment with non-vanishing interaction. We 
will refer to p^ as the thermal Green's function; this 
function can be expanded perturbatively in the system- 
environment coupling as 



Pp 



-/3H 



S P,3 



(74) 



Zoifi) 

where Za(j3) is the partition function of the free system 
The second-order corrections are given by 



(Ui\6pp \u)j) = 



l ijk 

.zW) 



Wk) (Wfe| CT Xn \Uj 



(75) 

All terms with u>i = loj are zero, so that this expression 
gives no correction to the diagonal elements of the den- 
sity matrix. Otherwise, the (non-resonant) off-diagonal 
coefficients are given by 



R?%\ =Im 



Afj 



i(wi/c)— A r , 

LU,, - Uj 



+ Im 



e~P Wi A mn (W« ) - e~^i Aran (u>kj ) 



LOi — UJj 



(76) 



with the free ground-state energy set to zero. These co- 
efficients agree perturbatively with those from Eq. (72). 
Because such an expansion is inherently secular in j3, it is 
valid only at a sufficiently high temperature such that the 
perturbations are small compared to the smallest Boltz- 
mann weight, 



n(0„,T) 



n(n n ,T) + lJ Vn(n m ,r) + i 



n(Q m ,T) 



7 „ -/3(n„+n ro ) 
^ 

The expansion does not apply at lower temperatures. Re- 
liability of the expansion at higher temperature suggests 
that the diagonal corrections to the asymptotic state 
must be suppressed there. 

Since neither the second-order master equation nor the 
perturbative expansion of the thermal Green's function 
can give the full low-temperature solution, including di- 
agonal corrections, it appears that in general this will re- 
quire the fourth-order master equation coefficients. How- 
ever, at zero temperature the thermal state is simply the 
ground state of the total system-environment Hamilto- 
nian. This ground state can be calculated perturbatively 
from the Hamiltonian as usual in a closed system, and the 
zero-temperature reduced thermal state follows directly. 
All three of these formalisms are fully consistent as shown 
in Ref. [17]. At zero temperature the off-diagonal second- 
order corrections to the asymptotic state are still of the 
form given in Eqs. (75) and (76), with the coefficients 
evaluated in the limit j3 — > oo. The diagonal (and reso- 
nant) perturbations are given by 



(78) 



lim Im 



■A n m{uik) + e 



-/3u>i 



did 



'Amn ij^ki) 



where only a handful of terms are non- vanishing. We 
note that the expression inside the limit in Eq. (78) has 
both the correct low and high-temperature limits, so it 
may be roughly correct for all temperatures, but we have 
yet to fully investigate the fourth-order master equation. 

For most regimes the second-order thermal state can 
now be expressed entirely in terms of the second-order 
master equation coefficients and limits thereof. Therefore 
we can say that the environmentally induced correlations 
do vanish for large separations with a power-law decay 
like 1/r and 1/r 2 . 



B. Entanglement of Two Atoms 

Now we will consider the bipartite entanglement be- 
tween any two atoms, labeled n and m in a common 
quantum field. We begin with some remarks that apply 
to any system of two qubits. We focus on the late-time 
dynamics of this system; we will compute the reduced 
density matrix for their asymptotic state p nm and derive 
the asymptotic value of entanglement between these two 
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atoms. We will see that this computation will also al- 
low us to show that all entangled initial states become 
disentangled at a finite time. 

To quantify the bipartite entanglement we will use 
Wootters' concurrence function [54], which is a mono- 
tone with a one-to-one relationship to the entanglement 
of formation for two qubits. The concurrence is defined 
as 

C(Pnm) =max{0,C(p nm )} (79) 
C( Pnm ) =^X~i -^~2-Vh-V^~4 (80) 

where Ai > A2 > A3 > A4 are the eigenvalues of the 
matrix 

Pnm Pnm = Pnm { a Vn a Vm Pnm a V n °2M> ) ) ( 81 ) 

which are always non-negative. A two-qubit state is en- 
tangled if and only if C > 0. It is important to note that 
C (p) is a continuous function of the matrix elements of p 
(since the eigenvalues of a matrix can be written as a con- 
tinuous function of the matrix elements [55]); this then 
implies that any density matrix with C_ < lies in the in- 
terior of the set of separable states (with every sufficiently 
nearby state also separable), while states with C > lie 
in the interior of the set of entangled states. States with 
C_ = are separable but include states that lie on the 
boundary between the two sets, infinitesimally close to 
both entangled states and the interior of the separable 
states. Any separable pure state lies on this boundary 
[56]. 

Given the late-time asymptotic state of two atoms 
p nm , one can easily compute the asymptotic entangle- 
ment from C_(p nm ). Based on the preceding paragraph, 
however, we know that this will also tell us something 
qualitatively about the late-time entanglement dynam- 
ics. If C (p nm ) < then (assuming only continuous evo- 
lution in state space) every initial state must become 
separable at some finite time as it crosses into the set of 
separable states. Likewise, if C(p nm ) > then all ini- 
tial states lead to entanglement at sufficiently late time 
and any sudden death of entanglement must be followed 
by revival. In models such as ours which have a unique 
asymptotic state, it is only when C(p nm ) = that this 
qualitative feature of the late-time behavior will depend 
on the initial state, with some entangled states remaining 
separable after some finite time and others becoming dis- 
tcntangled only asymptotically in the limit t — > 00 as in 
[7, 11]. Previous work has pointed out that the late-time 
entanglement dynamics can be determined by the asymp- 
totic state in this way [57, 58], with Yu and Eberly [59] 
discussing the role of C in predicting sudden death. In 
Refs. [58, 59] the authors consider models with multiple 
steady states, which introduces additional dependence on 
initial conditions. 

It can be seen that none of the foregoing discussion 
is specific to the concurrence; it would apply to any 
quantity that is a continuous function of the density ma- 
trix, takes on negative values for some separable states, 
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FIG. 8: Unmaximized concurrence for 2 atoms initially 
in the state t/> = \/0\T |0, 0) + \/CK9 |1, 1) with 
7 = J7/1000 and r — 2n/20il, showing sudden death, 
revival, and (seemingly) asymptotic disentanglement as 

the concurrence trends to zero. (Compare with [11] 
Fig. 2.) The middle, red region spans ±7/51 and denotes 
the order of neighborhood around C — which cannot 
be resolved by any second-order master equation. 

c 




FIG. 9: The same plot as Fig. 8, but with the larger 
decay rate 7 = fi/100. In this case the prediction of 
entanglement revival is swamped by uncertainties 
inherent in the second-order master equation, the order 
of which are spanned by the middle, red region. 



and is an entanglement monotone when non-negative. If 
we have such an unmaximized entanglement function £ 
from which an entanglement monotone can be defined 
by £ — max{0,£}, then we can use it just as we have 
discussed using C above. As illustrated with a specific 
example in Fig. 10, entanglement sudden death occurs 
because the unmaximized entanglement function asymp- 
totes towards a negative value, whereas any entanglement 
monotone (derived from £_ or otherwise) cannot go below 
zero, leading to an abrupt sudden death of entanglement 
when £_ becomes negative. 

An important point arises from the facts we have noted 
about C_ and separability: At sufficiently low temper- 
ature the 0(7) corrections to the asymptotic state are 
required to calculate the sign of C(p nm ) and, therefore, 
even the qualitative features of late-time entanglement 
dynamics. Specific examples are plotted in Fig. 8-9. At 
zero temperature, the zeroth-order asymptotic state is 
simply the ground state of the system (assuming no de- 
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FIG. 10: The same plot as Fig. 8, but at finite 
temperature with T = 0/5. In this case sudden death 
of entanglement (after a brief revival) can be confirmed 
with the second-order master equation, but only by 

using an unmaximized entanglement monotone. 



generacy at the ground energy) according to Eq. (70). 
So the zeroth-order asymptotic state is a pure separable 
state. This means that it lies on the boundary between 
the entangled and separable states, and in general some 
initial states will suffer sudden death while others will 
not, as depicted in Fig. 11a. But any non-zero pertur- 
bation, however small, can lead to asymptotic entangle- 
ment or can place the asymptotic state in the interior of 
the separable states, implying sudden death for all ini- 
tial conditions. Fig. lib shows each of these situations. 
Thus, knowing only the zeroth-order asymptotic state 
one can make no meaningful prediction about late-time 
entanglement dynamics, and this will always be the case 
when using the rotating-wave approximation, because it 
neglects the second-order corrections to the asymptotic 
state [ ]. This makes calculations such as [11] inca- 
pable of correctly predicting these features. Moreover, as 
one can see from Fig. 9, in certain regimes errors inher- 
ent in the second-order master equation can completely 
swamp almost all predictions of the entanglement dy- 
namics. For example, with the initial conditions that 
give rise to Fig. 8-9, one can see that the decay rate 
must satisfy 7 -C f2/100 in order to ensure that the pre- 
diction of revival, made for that situation in Fig. 2 of 
[11], is accurate despite the presence of errors. This is a 
general consequence of [1] that applies to any calculation 
resulting in a concurrence comparable to 7/O. 

At positive temperature the zeroth-order asymptotic 
state is simply the Boltzmann state p Tl which lies in the 
interior of the set of separable states [57], and 



PtPt 



-(fi„+n m )/r 



(82) 



so that C_{Pt) = -2e-(°™+ n -)/( 2T )/Z (T). The C( 7 ) 
corrections to p nm will yield order 0(7) corrections to 
PnmPnm- Then simply from the definition of C_ we know 
that so long as the temperature is sufficiently high that 
Eq. (77) is satisfied the corrections to p nm will cause at 
most 0(7) corrections to C (p nm ) so that it must remain 
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(a) Pure Asymptotic State 

Ent. 




(b) Mixed Asymptotic State 

FIG. 11: A schematic representation of the evolution in 
state space. The white area represents entangled states 
{C_> 0), while the gray areas represent separable states 

C_ < with the dark gray representing states with 
(7 = 0. The asymptotic state is represented by o, while 

initial states are represented by ■. In (a) we have the 
asymptotic state on the boundary as in the zeroth-order 
at T = 0. In (b) two scenarios are shown that can arise 
from a small perturbation moving the asymptotic state 
off the boundary, into the interior of one of the two sets. 

This illustrates how such a perturbation qualitatively 
changes the late-time entanglement dynamics. 



negative. Consequently, the second-order asymptotic 
state still lies in the interior of the separable states, and 
all initial states will suffer entanglement sudden death 
at sufficiently late times. For lower temperatures it does 
not appear that the sign of C(p nm ) can be generically 
predicted, and one must find the late-time asymptotic 
state for the specific system in question which generally 
requires terms from the fourth-order master equation. 
Returning to the specifics of the particular model ex- 
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amined in this paper, from Eq. (75) it can be readily 
seen that the atoms are correlated in the asymptotic 
state at all temperatures, and from our second-order co- 
efficients these correlations experience power-law decay 
with separation. However, we find based on Eqs. (74), 
(75), and (76) that when the high-temperature expan- 
sion is valid (according to Eq. (77)) the asymptotic state 
has Q(p nm ) < 0. At zero temperature, Eqs. (76) and 
(78) also give C_(p nm ) < 0. In both cases the asymp- 
totic state lies in the interior of the separable states, 
and all initial states become separable permanently af- 
ter some finite time. With this property upheld for zero 
and high temperatures, we suspect this to be true at 
all temperatures, making entanglement sudden death a 
generic feature which happens in every case in this model. 
Of course, as discussed in Sec. IV B, for closely spaced 
atoms there can be a dark state, so that entanglement 
persists over a long timescale before eventually succumb- 
ing to sudden death. It should also be noted that, while 
this examination of the asymptotic behavior tells us that 
entanglement always remains zero after some finite time, 
we do find 0(7°) sudden death and revival of entangle- 
ment at earlier times for some initial states (similar to 
[11]) as depicted in Fig. 8-10. 

In Fig. 12 we plot C as it varies with separation dis- 
tance and frequency detuning. As a consistency check 
we also calculated the logarithmic negativity and found 
it to be consistent with the concurrence to second order. 
The behavior of the entanglement is markedly different 
from that of two oscillators in a field. The separation 
dependence of two resonant oscillators was considered in 
Ref. [36] and the more general solution will be given in 
Ref. [60] . For two oscillators, there can be asymptotic en- 
tanglement if they are held very close and near enough to 
resonance with each other. Separation and detuning then 
causes the entanglement monotones to decay away. For 
the two-atom case studied here asymptotic entanglement 
does not exist, and resonant tuning with proximity will 
only exacerbate the problem. Permanent sudden death 
of entanglement occurs because the unmaximized entan- 
glement functions can trend below zero within a finite 
amount of time and without the need of any asymptotic 
limit. 



VI. DISCUSSION 

In this paper we have derived the dynamics of a collec- 
tion of two-level atoms under a dipole approximation in- 
teracting with a common quantized electromagnetic field 
assuming only weak coupling to the field. With a careful 
perturbative analysis we have obtained all dynamical in- 
formation to second order, including regimes which can- 
not be described by RWA-Lindblad equations. We have 
also presented a method of finding the zero-temperature 
asymptotic state to higher accuracy than is possible di- 
rectly with any second-order master equation, including 
those derived with the RWA and BMA. We have used 



this to show that even at zero temperature the bipar- 
tite entanglement between any pair of two-level atoms 
will undergo sudden death for all initial atomic states, in 
contrast to the predictions of previous theoretical treat- 
ments [11]. (We will point out specific deficiencies of 
[12] in a later communication.) Finally, we have noticed 
that a class of iV-atom (Dickie) dark states, which ap- 
pear normal in the RWA Hamiltonian, are not ordinary 
dark states and are, in fact, temperature dependent. 

We have argued that in the RWA there can be inac- 
curacies in all entries of the density matrix that are of 
the order of the damping rate 7. By contrast, when rep- 
resented in the (free) energy basis the solution we have 
derived here will have off-diagonal elements that are ac- 
curate at second-order, having 0(7 2 ) errors. Even in 
this solution diagonal matrix elements (and matrix el- 
ements between any two degenerate energy states) can 
still have 0(7) errors, due to a fundamental limitation 
of any second-order master equation. However, the ex- 
pectation of any operator that has vanishing diagonals 
in the energy basis (including atomic dipole operators), 
will have only 0(7 2 ) errors. Moreover, unlike some other 
methods of solution, our solution can be applied when 
the atoms have distinct (but close) frequencies. 

At sufficiently low temperature, the zeroth-order 
asymptotic state (given by the RWA) is near the bound- 
ary between the separable and entangled states, and 
the small perturbation induced by the environment at 
0(7) can push it into either set. Depending on which 
set the perturbed asymptotic states fall into, all states 
may experience entanglement sudden death or all may 
become entangled asymptotically. We have presented a 
second-order solution for the asymptotic state of any two 
atoms, which allows us to say decisively that the zero- 
temperature asymptotic state of those atoms is separable, 
and pairwise entanglement of all two-level atoms experi- 
ences sudden death regardless of the initial state. 

It should be noted that, for example, in some optical- 
frequency atomic systems the 0(7) corrections we discuss 
can be quite small, with "f/fl being perhaps something on 
the order of 10 -9 . Though lowest order corrections to the 
timescales cannot be ignored (as they are responsible for 
the presence of dissipation), corrections of this size to the 
values of the density matrix elements at any instant can 
easily be considered negligible. However, in the case of 
a theoretical study of entanglement sudden death, where 
one wishes to distinguish asymptotic decay to zero from 
vanishing in finite time, small perturbations can become 
vitally important, as they do at low temperature. And in 
optical frequency atomic systems at room temperature 
the thermal-average photon number will be far smaller 
than 10 -9 , placing the system deep into what we are 
considering the low-temperature regime for entanglement 
dynamics. 

We have characterized the sub- and super-radiant 
states that exist in this model when the RWA is not 
used. We have shown that there is still a long-lived, 
highly-entangled dark state when the atoms have small 
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FIG. 12: Unmaximized concurrence for two resonant atoms at various separation distances and frequency detunings 
at zero temperature and for 7 = (f2)/100. The left plot includes all terms in the model whereas the right plot 
neglects the 1/r magnetostatic interactions, as to inspect the finite contributions. For the typical RWA-Lindblad 
equation, these plots would simply be C_ = given asymptotic relaxation to the free ground state. The plots here 
are accurate to second order, which is beyond that of any prediction that can be made with a second-order master 

equation. 



enough separation, and sudden death of entanglement 
occurs only on the much longer timescale of decay of 
this state (assuming it had some population in the ini- 
tial state). In this simple model, decoherence-free dark 
states are achievable for arbitrary temperature and dis- 
sipation, whereas typically these factors together are the 
primary cause of decoherence. However, for N > 2, the 
finite-temperature Dickie dark states fall into two cat- 
egories: one group is temperature-independent and the 
other becomes temperature dependent. This tempera- 
ture dependence is not seen after employing the RWA. 

We close with a few remarks: 1) With the knowledge 
of distance dependence, to preserve entanglement in time 
one should place the atoms very close to each other in 
the field, so as to produce strong cross correlations in the 
noise. But at some proximity one must also consider fur- 
ther atom-atom interactions self-consistently within the 
confines of field theory. 2) Qualitative differences be- 



tween systems under the two-level and dipole-interaction 
approximations and harmonic systems suggests a degree 
of model dependence in some of the phenomena con- 
sidered; this merits further investigation into the con- 
sequences of these approximations. 3) Many other sorts 
of level structures are relevant to experimental systems, 
both in terms of the number of levels involved and the an- 
gular momentum exchange with the field. The method- 
ology and conceptions developed in this work can be ap- 
plied for the analysis of the non-Markovian dynamics of 
more general systems, from which one can perhaps bet- 
ter understand how model-dependent the entanglement 
behavior considered herewith is. 
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